********************************************************************************
********************************************************************************
*********************** RANKING ARCH. CONTRACTORS ACROSS PROGRAM YEARS *********

** load data
clear all
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"

** merge with info about furnace replacements
merge m:1 Household using "$maindir/Data/furn_replace.dta"
drop if _merge==2
drop _merge

***** some extra restrictions on number of homes served by each contractor
gen aux=1
sort Household meterend
bysort Household : gen homeobs = sum(aux) if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs = . if homeobs!=1

sort Household meterend
bysort Household : gen homeobs_post = sum(aux) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs_post = . if homeobs_post!=1

** fixing builddate to avoid omitted variables
replace builddate = 11 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

** fixing sqfeet to avoid omitted variables
drop roundsqfeet
gen roundsqfeet = .
	forval i = 0(500)2500 {
		replace roundsqfeet = `i' if sqfeet>`i' & sqfeet<=`i'+500
	}
replace roundsqfeet = 3000 if sqfeet>3000 & sqfeet!=.

** fixing some spending variables to avoid omitted bins
drop binAirCon
gen binAirCon = .
replace binAirCon=0 if Real_tot_actAirCon==0
replace binAirCon=1 if Real_tot_actAirCon>0 & Real_tot_actAirCon!=.

replace binAttic = 2700 if Real_tot_actAttic>2400 & Real_tot_actAttic!=.

** demographic variables
gen simpocc = noccupants
replace simpocc = 5 if noccupants>=5 & noccupants!=.

gen binAge = .
forval i = 20(10)80 {
replace binAge = `i' if Age>`i' & Age<=`i'+10
}

gen binIncome = .
forval i = 5000(5000)35000 {
replace binIncome = `i' if Real_income>`i' & Real_income<=`i'+5000
}
replace binIncome = 1 if Real_income<=5000
replace binIncome = 40000 if Real_income>40000 & Real_income!=.


*** housing structure
gen binSqft = .
forval i = 600(300)2700 {
replace binSqft = `i' if sqfeet>`i' & sqfeet<=`i'+300
}
replace binSqft = 1 if sqfeet<=600
replace binSqft = 3000 if sqfeet>3000 & sqfeet!=.

drop BlowerPreBins
gen BlowerPreBins = .
forval i = 2000(1000)9000 {
replace BlowerPreBins = `i' if Blower_Pre>`i' & Blower_Pre<=`i'+1000
}
replace BlowerPreBins = 1 if Blower_Pre<=2000
replace BlowerPreBins = 10000 if Blower_Pre>10000 & Blower_Pre!=.

*** energy prices
gen binelec = .
local counter = 2
forval i = 10.5(0.5)12.5 {
replace binelec = `counter' if elec_prices>`i' & elec_prices<=`i'+0.5
local counter = `counter' + 1
}
replace binelec = 1 if elec_prices<10.5
replace binelec = `counter' if elec_prices>13 & elec_prices!=.

gen bingas = .
local counter = 2
forval i = 7(1)16 {
replace bingas = `counter' if gas_prices>`i' & gas_prices<=`i'+0.5
local counter = `counter' + 1
}
replace bingas = 1 if gas_prices<7
replace bingas = `counter' if gas_prices>17 & gas_prices!=.


********************************************************************************
****** finally, regression for the wedge
* transforming decimals in percent
gen savings_gap_pct = savings_gap*100

forval i = 1/200 {
gen savings_gap_pct`i' = 100*savings_gap`i'
}

gen realized_savings_pct = realized_savings*100
winsor2 realized_savings_pct, trim replace cuts(1 99)

forval i = 1/200 {
gen realized_savings_pct`i' = 100*realized_savings`i'
winsor2 realized_savings_pct`i', trim replace cuts(1 99)
}


** main heat pre-treatment
gen MainHeatkBTU = MainHeatBTU/1000

gen binheatkbtu = .
forval i = 60(10)120 {
replace binheatkbtu = `i' if MainHeatkBTU>`i' & MainHeatkBTU<=`i'+10
}
replace binheatkbtu = 1 if MainHeatkBTU<=60
replace binheatkbtu = 130 if MainHeatkBTU>130 & MainHeatkBTU!=.

gen heatworks = 0 if MainHeatOperational=="False"
replace heatworks = 1 if MainHeatOperational=="True"


* fixing variables that have a zero indicator (for later interacting)
foreach x in nstories binAirCon white black hispanic otherrace haselderly hasminor heatworks {
replace `x' = `x'+1
}
replace roundAtticR = 1 if roundAtticR==0

gen heatnotwork = 2 if MainHeatOperational=="False"
replace heatnotwork = 1 if MainHeatOperational=="True"

replace binFurnace = 1 if binFurnace==0

***** some missing vars
gen pct_BlowerReduc = BlowerReduc/Blower_Pre

gen atticR = Attic_RValue
replace atticR = . if atticR==0

gen hasattic = 0 if Attic_RValue==0
replace hasattic = 1 if Attic_RValue>0 & Attic_RValue!=.

winsor2 sqfeet, replace trim cuts(1 99)

gen BlowerReduc100 = BlowerReduc/100

gen summer = 0
replace summer = 1 if end_month>=4 & end_month<=9

***** climate zones
merge m:1 Household using "$maindir/Data/wx_weather.dta", keepusing(region)
drop if _merge==2
drop _merge
egen regionID = group(region)

*** sample to be used in regressions
reghdfe log_total_mmbtu i.treated c.HDD60 c.HDD60#c.HDD60 c.CDD75 c.CDD75#c.CDD75 c.precip c.tmax c.tmin if twoyears_prepost==1 & meterdays_monthly!=., absorb(i.Household i.month_year) vce(cluster Household end_month)
drop regsample
gen regsample = e(sample) /* to keep stable sample across regressions */
gen postsample = regsample if treated==1

* observations post
drop homeobs_post
sort Household postsample
bysort Household : gen homeobs_post = sum(aux) if postsample==1
replace homeobs_post = . if homeobs_post!=1

* observations by agency
sort AgencyID postsample
bysort AgencyID : gen ageobs = sum(aux) if postsample==1
bysort AgencyID : egen Nagency = max(ageobs)
drop ageobs

merge m:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(arch_contID) nogen keep(3)
* contractors by agency
forval i = 2009/2016 {
sort AgencyID arch_contID postsample
bysort AgencyID arch_contID : gen firstobs = sum(aux) if postsample==1 & ProgramYear==`i'
replace firstobs = . if firstobs!=1
bysort AgencyID : gen maxobs = sum(firstobs)
bysort AgencyID : egen agency_cont_obs`i' = max(maxobs)
replace agency_cont_obs`i' = . if AgencyID==.
drop firstobs maxobs
}

* jobs by contractor
sort arch_contID Household postsample
bysort arch_contID Household : gen firstobs = sum(aux) if postsample==1
replace firstobs = . if firstobs!=1
bysort arch_contID : gen maxobs = sum(firstobs)
bysort arch_contID : egen cont_jobs = max(maxobs)
replace cont_jobs = . if arch_contID==.
drop firstobs maxobs

* jobs by contractor and by PY
sort arch_contID Household postsample

forval i = 2009/2016 {
bysort arch_contID Household : gen firstobs = sum(aux) if postsample==1 & ProgramYear==`i'
replace firstobs = . if firstobs!=1
bysort arch_contID : gen maxobs = sum(firstobs)
bysort arch_contID : egen cont_jobs_py`i' = max(maxobs)
replace cont_jobs_py`i' = . if arch_contID==.
drop firstobs maxobs
}

forval i = 2009/2016 {
tab arch_contID if postsample==1 & ProgramYear==`i' & treated==1 & cont_jobs_py`i'>=10 & cont_jobs_py`i'!=. & regionID!=1
di `r(r)'
}

****** yearly contractor fixed effects
** specification with architectural contractors only
set matsize 10000


** new definition for bins of Wall Insulation
replace binWallIns = 2400 if binWallIns>=2400 & binWallIns!=.


**** alternative (coarse) binned spending for interactions
gen coarsebinAirSeal = 1 if binAirSeal==1
replace coarsebinAirSeal = 2 if binAirSeal>1 & binAirSeal<=300
replace coarsebinAirSeal = 3 if binAirSeal>300 & binAirSeal<=600
replace coarsebinAirSeal = 4 if binAirSeal>600 & binAirSeal!=.

gen coarsebinWallIns = 1 if binWallIns==1
replace coarsebinWallIns = 2 if binWallIns>1 & binWallIns<=600
replace coarsebinWallIns = 3 if binWallIns>600 & binWallIns<=1200
replace coarsebinWallIns = 4 if binWallIns>1200 & binWallIns<=1800
replace coarsebinWallIns = 5 if binWallIns>1800 & binWallIns<=2400
replace coarsebinWallIns = 6 if binWallIns>2400 & binWallIns!=.

gen coarsebinAttic = 1 if binAttic==1
replace coarsebinAttic = 2 if binAttic>1 & binAttic<=600
replace coarsebinAttic = 3 if binAttic>600 & binAttic<=1200
replace coarsebinAttic = 4 if binAttic>1200 & binAttic<=1800
replace coarsebinAttic = 5 if binAttic>1800 & binAttic<=2400
replace coarsebinAttic = 6 if binAttic>2400 & binAttic!=.

gen coarsebinFoundation = 1 if binFoundation==1
replace coarsebinFoundation = 2 if binFoundation>1 & binFoundation<=400
replace coarsebinFoundation = 3 if binFoundation>400 & binFoundation<=800
replace coarsebinFoundation = 4 if binFoundation>800 & binFoundation<=1200
replace coarsebinFoundation = 5 if binFoundation>1200 & binFoundation<=1600
replace coarsebinFoundation = 6 if binFoundation>1600 & binFoundation!=.

gen coarsebinWindow = 1 if binWindow==1
replace coarsebinWindow = 2 if binWindow>1 & binWindow<=400
replace coarsebinWindow = 3 if binWindow>400 & binWindow<=800
replace coarsebinWindow = 4 if binWindow>800 & binWindow<=1200
replace coarsebinWindow = 5 if binWindow>1200 & binWindow<=1600
replace coarsebinWindow = 6 if binWindow>1600 & binWindow!=.




********************************************************************************
******************************** regressions for contractor quality
reg realized_savings_pct ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  ///
		ib2009.ProgramYear##ib236.arch_contID ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID ib3000.BlowerPreBins ib1.builddate ///
		i.nstories ib1500.binSqft ///
		if regsample==1 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. & ProgramYear>=2009
parmest, saving("$maindir/Results/Model_Outputs/contfe_save_pooled.dta", replace) omit empty


reg realized_savings_pct ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib1.coarsebinAirSeal##ib1.coarsebinWallIns ///
		ib1.coarsebinAirSeal##ib1.coarsebinAttic ///
		ib1.coarsebinAirSeal##ib1.coarsebinFoundation ///	
		ib1.coarsebinAirSeal##ib1.coarsebinWindow ///	
		ib1.coarsebinAirSeal##ib3000.BlowerPreBins ///	
		ib1.coarsebinWallIns##ib1.coarsebinAttic ///
		ib1.coarsebinWallIns##ib1.coarsebinFoundation ///	
		ib1.coarsebinWallIns##ib1.coarsebinWindow ///	
		ib1.coarsebinWallIns##ib3000.BlowerPreBins ///	
		ib1.coarsebinAttic##ib1.coarsebinFoundation ///	
		ib1.coarsebinAttic##ib1.coarsebinWindow ///	
		ib1.coarsebinAttic##ib3000.BlowerPreBins ///	
		ib1.coarsebinFoundation##ib1.coarsebinWindow ///	
		ib1.coarsebinFoundation##ib3000.BlowerPreBins ///	
		ib1.coarsebinWindow##ib3000.BlowerPreBins ///	
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  ///
		ib2009.ProgramYear##ib236.arch_contID ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID##ib3000.BlowerPreBins i.regionID##ib1.builddate ///
		i.regionID##i.nstories i.regionID##ib1500.binSqft ///
		if regsample==1 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. & ProgramYear>=2009
gen qualregsamp_comp = e(sample)
parmest, saving("$maindir/Results/Model_Outputs/contfe_save_complex.dta", replace) omit empty



******************************
******** merge with contractor ranking by savings - pooled contractor regression - quality based on savings
tempfile ihwapgap
save `ihwapgap'
		
clear
use "$maindir/Results/Model_Outputs/contfe_save_pooled.dta"

forval i = 10/16 {
	gen py`i' = estimate if parm=="20`i'.ProgramYear"
	egen py`i'est = max(py`i')
	drop py`i'
}

keep if substr(parm, -5, .)=="ontID"

gen arch_contID = subinstr(parm, ".arch_contID", "", .)
forval i = 2009/2016 {
	replace arch_contID = subinstr(arch_contID, "`i'b.ProgramYear#", "", .)
	replace arch_contID = subinstr(arch_contID, "`i'o.ProgramYear#", "", .)
	replace arch_contID = subinstr(arch_contID, "`i'.ProgramYear#", "", .)
}

replace arch_contID = subinstr(arch_contID, "b", "", .)
replace arch_contID = subinstr(arch_contID, "o", "", .)
destring arch_contID, replace
	
gen ProgramYear = substr(par, 1, 4)
replace ProgramYear = subinstr(ProgramYear, ".", "", .)
replace ProgramYear = subinstr(ProgramYear, "a", "", .)
replace ProgramYear = subinstr(ProgramYear, "b", "", .)
replace ProgramYear = subinstr(ProgramYear, "o", "", .)
replace ProgramYear = subinstr(ProgramYear, "r", "", .)
destring ProgramYear, replace
replace ProgramYear = 2009 if ProgramYear<=2009
drop if substr(parm, 1, 5)=="2009b"

gen newest = estimate
replace omit = 0 if ProgramYear==2009 & arch_contID==236
forval i = 10/16 {
	replace omit = 0 if ProgramYear==20`i' & arch_contID==236
	replace newest = py`i'est if ProgramYear==20`i' & arch_contID==236
}
	
drop if omit==1 | empty==1
drop if newest==.

sort arch_contID ProgramYear

keep arch_contID ProgramYear newest
duplicates drop

rename newest contquality

xtset arch_contID ProgramYear
gen lag_contquality = L.contquality

bysort arch_contID : replace lag_contquality = contquality[_n-1] if lag_contquality==. & contquality[_n-1]!=.

*** ranking in percentiles instead of level variables
xtile contquality_rank = contquality  , n(5)

merge 1:m arch_contID ProgramYear using `ihwapgap'
drop _merge



******************************
******** merge with contractor ranking by savings - pooled contractor regression - quality based on savings - complex interactions
tempfile ihwapgap
save `ihwapgap'
		
clear
use "$maindir/Results/Model_Outputs/contfe_save_complex.dta"

forval i = 10/16 {
	gen py`i' = estimate if parm=="20`i'.ProgramYear"
	egen py`i'est = max(py`i')
	drop py`i'
}

keep if substr(parm, -5, .)=="ontID"

gen arch_contID = subinstr(parm, ".arch_contID", "", .)
forval i = 2009/2016 {
	replace arch_contID = subinstr(arch_contID, "`i'b.ProgramYear#", "", .)
	replace arch_contID = subinstr(arch_contID, "`i'o.ProgramYear#", "", .)
	replace arch_contID = subinstr(arch_contID, "`i'.ProgramYear#", "", .)
}

replace arch_contID = subinstr(arch_contID, "b", "", .)
replace arch_contID = subinstr(arch_contID, "o", "", .)
destring arch_contID, replace
	
gen ProgramYear = substr(par, 1, 4)
replace ProgramYear = subinstr(ProgramYear, ".", "", .)
replace ProgramYear = subinstr(ProgramYear, "a", "", .)
replace ProgramYear = subinstr(ProgramYear, "b", "", .)
replace ProgramYear = subinstr(ProgramYear, "o", "", .)
replace ProgramYear = subinstr(ProgramYear, "r", "", .)
destring ProgramYear, replace
replace ProgramYear = 2009 if ProgramYear<=2009
drop if substr(parm, 1, 5)=="2009b"

gen newest = estimate
replace omit = 0 if ProgramYear==2009 & arch_contID==236
forval i = 10/16 {
	replace omit = 0 if ProgramYear==20`i' & arch_contID==236
	replace newest = py`i'est if ProgramYear==20`i' & arch_contID==236
}
	
drop if omit==1 | empty==1
drop if newest==.

sort arch_contID ProgramYear

keep arch_contID ProgramYear newest
duplicates drop

rename newest contquality_comp

xtset arch_contID ProgramYear
gen lag_contquality_comp = L.contquality_comp

bysort arch_contID : replace lag_contquality_comp = contquality_comp[_n-1] if lag_contquality_comp==. & contquality_comp[_n-1]!=.

*** ranking in percentiles instead of level variables
xtile contquality_comp_rank = contquality_comp  , n(5)

merge 1:m arch_contID ProgramYear using `ihwapgap'
drop _merge



********************************************************************************
***************** Bootstrapping contractor quality measure

**** bootstrapping regressions for quality FE
forval i = 1/200 {

	di "starting bootstrap iteration `i'"
	
	qui : reg realized_savings_pct`i' ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
			ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
			ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
			ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1.roundAtticR  ///
			ib2009.ProgramYear##ib236.arch_contID ///
			c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
			i.end_month elec_prices gas_prices ///
			i.regionID ib3000.BlowerPreBins ib1.builddate ///
			i.nstories ib1500.binSqft [fweight = weight_boots`i'] ///
			if regsample==1 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. & ProgramYear>=2009
	parmest, saving("$maindir/Results/Model_Outputs/contfe_save_pooled`i'.dta", replace) omit empty

	qui : reg realized_savings_pct`i' ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
			ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
			ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
			ib1.coarsebinAirSeal##ib1.coarsebinWallIns ///
			ib1.coarsebinAirSeal##ib1.coarsebinAttic ///
			ib1.coarsebinAirSeal##ib1.coarsebinFoundation ///	
			ib1.coarsebinAirSeal##ib1.coarsebinWindow ///	
			ib1.coarsebinAirSeal##ib3000.BlowerPreBins ///	
			ib1.coarsebinWallIns##ib1.coarsebinAttic ///
			ib1.coarsebinWallIns##ib1.coarsebinFoundation ///	
			ib1.coarsebinWallIns##ib1.coarsebinWindow ///	
			ib1.coarsebinWallIns##ib3000.BlowerPreBins ///	
			ib1.coarsebinAttic##ib1.coarsebinFoundation ///	
			ib1.coarsebinAttic##ib1.coarsebinWindow ///	
			ib1.coarsebinAttic##ib3000.BlowerPreBins ///	
			ib1.coarsebinFoundation##ib1.coarsebinWindow ///	
			ib1.coarsebinFoundation##ib3000.BlowerPreBins ///	
			ib1.coarsebinWindow##ib3000.BlowerPreBins ///	
			ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1.roundAtticR  ///
			ib2009.ProgramYear##ib236.arch_contID ///
			c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
			i.end_month elec_prices gas_prices ///
			i.regionID##ib3000.BlowerPreBins i.regionID##ib1.builddate ///
			i.regionID##i.nstories i.regionID##ib1500.binSqft [fweight = weight_boots`i'] ///
			if regsample==1 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. & ProgramYear>=2009
	parmest, saving("$maindir/Results/Model_Outputs/contfe_save_complex`i'.dta", replace) omit empty

}



******************************
******** merge with contractor ranking by savings - pooled contractor regression - quality based on savings
tempfile ihwapgap
save `ihwapgap'
keep arch_contID ProgramYear
duplicates drop
tempfile contcoefs
save `contcoefs'

forval j = 1/200 {
		
	di "starting iteration `j'"
	
	clear
	use "$maindir/Results/Model_Outputs/contfe_save_pooled`j'.dta"

	forval i = 10/16 {
		gen py`i' = estimate if parm=="20`i'.ProgramYear"
		egen py`i'est = max(py`i')
		drop py`i'
	}

	keep if substr(parm, -5, .)=="ontID"

	gen arch_contID = subinstr(parm, ".arch_contID", "", .)
	forval i = 2009/2016 {
		replace arch_contID = subinstr(arch_contID, "`i'b.ProgramYear#", "", .)
		replace arch_contID = subinstr(arch_contID, "`i'o.ProgramYear#", "", .)
		replace arch_contID = subinstr(arch_contID, "`i'.ProgramYear#", "", .)
	}

	replace arch_contID = subinstr(arch_contID, "b", "", .)
	replace arch_contID = subinstr(arch_contID, "o", "", .)
	destring arch_contID, replace
		
	gen ProgramYear = substr(par, 1, 4)
	replace ProgramYear = subinstr(ProgramYear, ".", "", .)
	replace ProgramYear = subinstr(ProgramYear, "a", "", .)
	replace ProgramYear = subinstr(ProgramYear, "b", "", .)
	replace ProgramYear = subinstr(ProgramYear, "o", "", .)
	replace ProgramYear = subinstr(ProgramYear, "r", "", .)
	destring ProgramYear, replace
	replace ProgramYear = 2009 if ProgramYear<=2009
	drop if substr(parm, 1, 5)=="2009b"

	gen newest = estimate
	replace omit = 0 if ProgramYear==2009 & arch_contID==236
	forval i = 10/16 {
		replace omit = 0 if ProgramYear==20`i' & arch_contID==236
		replace newest = py`i'est if ProgramYear==20`i' & arch_contID==236
	}
		
	drop if omit==1 | empty==1
	drop if newest==.

	sort arch_contID ProgramYear

	keep arch_contID ProgramYear newest
	duplicates drop

	rename newest contquality

	xtset arch_contID ProgramYear
	gen lag_contquality = L.contquality

	bysort arch_contID : replace lag_contquality = contquality[_n-1] if lag_contquality==. & contquality[_n-1]!=.
	
	keep arch_contID ProgramYear contquality lag_contquality
	rename contquality contquality`j'
	rename lag_contquality lag_contquality`j'

	merge 1:1 arch_contID ProgramYear using `contcoefs', nogen
	save `contcoefs', replace

	******************************
	******** merge with contractor ranking by savings - pooled contractor regression - quality based on savings - complex interactions
		
	clear
	use "$maindir/Results/Model_Outputs/contfe_save_complex`j'.dta"

	forval i = 10/16 {
		gen py`i' = estimate if parm=="20`i'.ProgramYear"
		egen py`i'est = max(py`i')
		drop py`i'
	}

	keep if substr(parm, -5, .)=="ontID"

	gen arch_contID = subinstr(parm, ".arch_contID", "", .)
	forval i = 2009/2016 {
		replace arch_contID = subinstr(arch_contID, "`i'b.ProgramYear#", "", .)
		replace arch_contID = subinstr(arch_contID, "`i'o.ProgramYear#", "", .)
		replace arch_contID = subinstr(arch_contID, "`i'.ProgramYear#", "", .)
	}

	replace arch_contID = subinstr(arch_contID, "b", "", .)
	replace arch_contID = subinstr(arch_contID, "o", "", .)
	destring arch_contID, replace
		
	gen ProgramYear = substr(par, 1, 4)
	replace ProgramYear = subinstr(ProgramYear, ".", "", .)
	replace ProgramYear = subinstr(ProgramYear, "a", "", .)
	replace ProgramYear = subinstr(ProgramYear, "b", "", .)
	replace ProgramYear = subinstr(ProgramYear, "o", "", .)
	replace ProgramYear = subinstr(ProgramYear, "r", "", .)
	destring ProgramYear, replace
	replace ProgramYear = 2009 if ProgramYear<=2009
	drop if substr(parm, 1, 5)=="2009b"

	gen newest = estimate
	replace omit = 0 if ProgramYear==2009 & arch_contID==236
	forval i = 10/16 {
		replace omit = 0 if ProgramYear==20`i' & arch_contID==236
		replace newest = py`i'est if ProgramYear==20`i' & arch_contID==236
	}
		
	drop if omit==1 | empty==1
	drop if newest==.

	sort arch_contID ProgramYear

	keep arch_contID ProgramYear newest
	duplicates drop

	rename newest contquality_comp

	xtset arch_contID ProgramYear
	gen lag_contquality_comp = L.contquality_comp

	bysort arch_contID : replace lag_contquality_comp = contquality_comp[_n-1] if lag_contquality_comp==. & contquality_comp[_n-1]!=.
	
	keep arch_contID ProgramYear contquality_comp lag_contquality_comp
	rename contquality_comp contquality_comp`j'
	rename lag_contquality_comp lag_contquality_comp`j'

	merge 1:1 arch_contID ProgramYear using `contcoefs', nogen
	save `contcoefs', replace

	di "end of iteration `j'"
}
merge 1:m arch_contID ProgramYear using `ihwapgap'
drop _merge



********************************************************************************
********************************************************************************
******* new measure of quality

*** regress quality on observables plus lagged quality
reg contquality lag_contquality ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID ib3000.BlowerPreBins ib1.builddate ///
		i.nstories ib1500.binSqft ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & ProgramYear<=2015 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. ///		
		, cluster(Household) 	
	
est sto qualreg
matrix contcoefs_simp = e(b)
sca rsq_simp = e(r2)
sca nobs = e(N)
gen contsamp = e(sample)
matrix qualcoefs = e(b)
predict contqual_pred_new

sort Household contsamp
bysort Household : gen conthomes = sum(aux) if contsamp==1
replace conthomes = . if conthomes!=1

		
*** complex interactions
reg contquality_comp lag_contquality_comp ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib1.coarsebinAirSeal##ib1.coarsebinWallIns ///
		ib1.coarsebinAirSeal##ib1.coarsebinAttic ///
		ib1.coarsebinAirSeal##ib1.coarsebinFoundation ///	
		ib1.coarsebinAirSeal##ib1.coarsebinWindow ///	
		ib1.coarsebinAirSeal##ib3000.BlowerPreBins ///	
		ib1.coarsebinWallIns##ib1.coarsebinAttic ///
		ib1.coarsebinWallIns##ib1.coarsebinFoundation ///	
		ib1.coarsebinWallIns##ib1.coarsebinWindow ///	
		ib1.coarsebinWallIns##ib3000.BlowerPreBins ///	
		ib1.coarsebinAttic##ib1.coarsebinFoundation ///	
		ib1.coarsebinAttic##ib1.coarsebinWindow ///	
		ib1.coarsebinAttic##ib3000.BlowerPreBins ///	
		ib1.coarsebinFoundation##ib1.coarsebinWindow ///	
		ib1.coarsebinFoundation##ib3000.BlowerPreBins ///	
		ib1.coarsebinWindow##ib3000.BlowerPreBins ///	
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.end_month i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		elec_prices gas_prices ///
		i.regionID##ib3000.BlowerPreBins i.regionID##ib1.builddate ///
		i.regionID##i.nstories i.regionID##ib1500.binSqft ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & ProgramYear<=2015 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. ///
		, cluster(Household) 		
		
est sto qualreg_comp
matrix contcoefs_comp = e(b)
sca rsq_comp = e(r2)
gen contsamp_comp = e(sample)
matrix qualcoefs_comp = e(b)
predict contqual_pred_comp


************** bootstrapping first-stage regresions
*** regress quality on observables plus lagged quality

matrix bootscoefs_simp = (.)
matrix bootscoefs_comp = (.)

forval i = 1/200 {
	
	di "starting iteration `i'"

	qui : reg contquality`i' lag_contquality`i' ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
			ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
			ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
			ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1.roundAtticR  i.ProgramYear ///
			c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
			i.end_month elec_prices gas_prices ///
			i.regionID ib3000.BlowerPreBins ib1.builddate ///
			i.nstories ib1500.binSqft ///
			[fweight = weight_boots`i'] if contsamp_comp==1 ///		
			, cluster(Household) 	
		
	gen bootsample = e(sample)	
	matrix tempmat = e(b)
	matrix tempcoef = tempmat[1,1]
	matrix bootscoefs_simp = (bootscoefs_simp \ tempcoef)
	predict contqualpred_boots`i' if bootsample==1
			
	*** complex interactions
	qui : reg contquality_comp`i' lag_contquality_comp`i' ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
			ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
			ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
			ib1.coarsebinAirSeal##ib1.coarsebinWallIns ///
			ib1.coarsebinAirSeal##ib1.coarsebinAttic ///
			ib1.coarsebinAirSeal##ib1.coarsebinFoundation ///	
			ib1.coarsebinAirSeal##ib1.coarsebinWindow ///	
			ib1.coarsebinAirSeal##ib3000.BlowerPreBins ///	
			ib1.coarsebinWallIns##ib1.coarsebinAttic ///
			ib1.coarsebinWallIns##ib1.coarsebinFoundation ///	
			ib1.coarsebinWallIns##ib1.coarsebinWindow ///	
			ib1.coarsebinWallIns##ib3000.BlowerPreBins ///	
			ib1.coarsebinAttic##ib1.coarsebinFoundation ///	
			ib1.coarsebinAttic##ib1.coarsebinWindow ///	
			ib1.coarsebinAttic##ib3000.BlowerPreBins ///	
			ib1.coarsebinFoundation##ib1.coarsebinWindow ///	
			ib1.coarsebinFoundation##ib3000.BlowerPreBins ///	
			ib1.coarsebinWindow##ib3000.BlowerPreBins ///	
			ib70.binheatkbtu i.heatworks ///
			ib2.simpocc ib50.binAge ib15000.binIncome ///
			i.white i.black i.haselderly i.hasminor ///
			ib1.roundAtticR  i.end_month i.ProgramYear ///
			c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
			elec_prices gas_prices ///
			i.regionID##ib3000.BlowerPreBins i.regionID##ib1.builddate ///
			i.regionID##i.nstories i.regionID##ib1500.binSqft ///
			[fweight = weight_boots`i'] if contsamp_comp==1  ///
			, cluster(Household) 		
			
	matrix tempmat = e(b)
	matrix tempcoef = tempmat[1,1]
	matrix bootscoefs_comp = (bootscoefs_comp \ tempcoef)
	predict contqualpred_comp_boots`i' if bootsample==1
	drop bootsample

}


matrix bootscoefs_simp = bootscoefs_simp[2..., 1]
matrix bootscoefs_comp = bootscoefs_comp[2..., 1]

mata : st_matrix("cont_var_simp", mm_colvar(st_matrix("bootscoefs_simp")))
matrix cont_var_simp = diag(cont_var_simp)
mata : st_matrix("cont_var_comp", mm_colvar(st_matrix("bootscoefs_comp")))
matrix cont_var_comp = diag(cont_var_comp)

matrix conts_simp = contcoefs_simp[1,1]
matrix conts_comp = contcoefs_comp[1,1]
mat colnames conts_simp="lag_quality"
mat colnames conts_comp="lag_quality"
mat colnames cont_var_simp="lag_quality"
mat rownames cont_var_simp="lag_quality"
mat colnames cont_var_comp="lag_quality"
mat rownames cont_var_comp="lag_quality"

do "$maindir/Code/bootstrapres.do"
bootstrapres conts_simp cont_var_simp nobs
estadd scalar r2 = rsq_simp
est sto qualreg

bootstrapres conts_comp cont_var_comp nobs
estadd scalar r2 = rsq_comp
est sto qualreg_comp

esttab  qualreg qualreg_comp ///
				using "$maindir/Results/Tables/contqual_1ststage_save.tex", ///
				replace b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
				nonotes nonumbers rename(lag_quality a) ///	
				keep(a) coeflabels(a "Lagged Quality") ///
				stats(N r2, fmt(%9.0gc %9.4f) ///
				labels("Observations" "R-squared"))


********************************************************************************
****************** Oster test

****** fully saturated model with complex interactions
reg contquality lag_contquality ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib1.coarsebinAirSeal##ib1.coarsebinWallIns ///
		ib1.coarsebinAirSeal##ib1.coarsebinAttic ///
		ib1.coarsebinAirSeal##ib1.coarsebinFoundation ///	
		ib1.coarsebinAirSeal##ib1.coarsebinWindow ///	
		ib1.coarsebinAirSeal##ib3000.BlowerPreBins ///	
		ib1.coarsebinWallIns##ib1.coarsebinAttic ///
		ib1.coarsebinWallIns##ib1.coarsebinFoundation ///	
		ib1.coarsebinWallIns##ib1.coarsebinWindow ///	
		ib1.coarsebinWallIns##ib3000.BlowerPreBins ///	
		ib1.coarsebinAttic##ib1.coarsebinFoundation ///	
		ib1.coarsebinAttic##ib1.coarsebinWindow ///	
		ib1.coarsebinAttic##ib3000.BlowerPreBins ///	
		ib1.coarsebinFoundation##ib1.coarsebinWindow ///	
		ib1.coarsebinFoundation##ib3000.BlowerPreBins ///	
		ib1.coarsebinWindow##ib3000.BlowerPreBins ///	
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.end_month i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		elec_prices gas_prices ///
		i.regionID##ib3000.BlowerPreBins i.regionID##ib1.builddate ///
		i.regionID##i.nstories i.regionID##ib1500.binSqft ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & ProgramYear<=2015 & treated==1 & regionID!=1 & Nagency>1500 & Nagency!=. ///
		, cluster(Household) 
		
est sto satmodel
		
sca rsq1 = round(`e(r2)'*1.3,.0001)
sca rsq2 = round(`e(r2)'*1.5,.0001)
sca rsq3 = round(`e(r2)'*1.7,.0001)


*** simpler model
local allways_controls ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID ib3000.BlowerPreBins ib1.builddate ///
		i.nstories ib1500.binSqft


matrix results = (., .)		
foreach x in rsq1 rsq2 rsq3 {

	est restore satmodel

	psacalc  beta lag_contquality , mcontrol(`allways_controls') rmax(`=scalar(`x')')

	sca beta_corr = r(beta)
	sca diff = r(dist)
	matrix tempres = (beta_corr, diff)
	
	matrix results = (results \ tempres)
	
}		


matrix results_delta = (.)		
foreach x in rsq1 rsq2 rsq3 {

	est restore satmodel

	psacalc  delta lag_contquality , mcontrol(`allways_controls') rmax(`=scalar(`x')')

	sca delta = r(delta)

	matrix results_delta = (results_delta \ delta)
	
}		

matrix results = results[2..., 1...]
matrix results_delta = results_delta[2..., 1...]

mat colnames results = "Beta" "Sq Diff"
mat colnames results_delta = "Delta"

mat rownames results="R-Sq Max = `=round(scalar(rsq1), 0.01)'" "R-Sq Max = `=round(scalar(rsq2), 0.01)'" "R-Sq Max = `=round(scalar(rsq3), 0.01)'" 
mat rownames results_delta="R-Sq Max = `=round(scalar(rsq1), 0.01)'" "R-Sq Max = `=round(scalar(rsq2), 0.01)'" "R-Sq Max = `=round(scalar(rsq3), 0.01)'" 

matrix list results
matrix list results_delta
matrix results_comb = (results[1..., 1], results_delta)
matrix list results_comb

estout matrix(results_comb, fmt(%9.4fc)) using "$maindir/Results/Tables/oster_combined_results.tex", replace style(tex)



********************************************************************************
**** ranking contractors based on predicted quality

sort arch_contID ProgramYear
bysort arch_contID ProgramYear : egen mean_contqual = mean(contqual_pred_new)
order mean_contqual, after(contqual_pred_new)

bysort arch_contID ProgramYear : gen tagcont = 1 if _n==1

** different ranking by PY
forval i = 2010/2015 {
	xtile contqual_quint_py`i' = mean_contqual if contsamp==1 & ProgramYear==`i', n(5)
}

gen new_contqual_quintile = contqual_quint_py2010 if ProgramYear==2010
forval i = 2011/2015 {
	replace new_contqual_quintile = contqual_quint_py`i' if ProgramYear==`i'
}
rename new_contqual_quintile contqual_quintile



********************************************************************************
********* HISTOGRAMS FOR DOLLARS SPENT IN EACH BIN / CONTRACTOR QUALITY

**** histograms of spending for contractor sample
label var Real_tot_actAirCon "Air Conditioning"
label var Real_tot_actAirSeal "Air Sealing"
label var Real_tot_actAttic "Attic"
label var Real_tot_actBaseload "Baseload"
label var Real_tot_actDoor "Door"
label var Real_tot_actFoundation "Foundation"
label var Real_tot_actFurnace "Furnace"
label var Real_tot_actGeneral "General"
label var Real_tot_actHealSfty "Health and Safety"
label var Real_tot_actWallIns "Wall Insulation"
label var Real_tot_actWindow "Window"
label var Real_tot_actWtHtr "Water Heater"

gen newbinAirSeal = binAirSeal - 100
replace newbinAirSeal = -100 if newbinAirSeal==-99
gen newbinAttic = binAttic - 300
replace newbinAttic = -300 if newbinAttic==-299
gen newbinWallIns = binWallIns - 300
replace newbinWallIns = -300 if newbinWallIns==-299
gen newbinFurnace = binFurnace - 300
replace newbinFurnace = -300 if newbinFurnace==-299
gen newbinWindow = binWindow - 200
replace newbinWindow = -200 if newbinWindow==-199

label define furnacelab2 -300 "0" 0 "1-300" 300 "301-600" 600 "601-900" 900 "901-1200" 1200 "1201-1500" ///
		1500 "1501-1800" 1800 "1801-2100" 2100 "2101-2400" 2400 "2401-2700" 2700 "2701-3000" 3000 ">3000"
label values newbinFurnace furnacelab2

label define wallinslab2 -300 "0" 0 "1-300" 300 "301-600" 600 "601-900" 900 "901-1200" 1200 "1201-1500" ///
		1500 "1501-1800" 1800 "1801-2100" 2100 ">2100"
label values newbinWallIns wallinslab2

label define airseallab2 -100 "0" 0 "1-100" 100 "101-200" 200 "201-300" 300 "301-400" 400 "401-500" ///
		500 "501-600" 600 "601-700" 700 "701-800" 800 "801-900" 900 ">900"
label values newbinAirSeal airseallab2
		
label define atticlab2 -300 "0" 0 "1-300" 300 "301-600" 600 "601-900" 900 "901-1200" 1200 "1201-1500" ///
		1500 "1501-1800" 1800 "1801-2100" 2100 "2101-2400" 2400 ">2400"
label values newbinAttic atticlab2

label define windowlab2 -200 "0" 0 "1-200" 200 "201-400" 400 "401-600" 600 "601-800" 800 "801-1000" ///
		1000 "1001-1200" 1200 "1201-1400" 1400 "1401-1600" 1600 "1601-1800" 1800 "1801-2000" 2000 ">2000"
label values newbinWindow windowlab2


foreach x of varlist binFurnace binAirSeal binAttic binWallIns binWindow {
	bysort new`x' conthomes : gen count`x' = _n if conthomes==1
	bysort new`x' : egen freq_`x'_conthomes = max(count`x') if conthomes==1
	drop count`x'
}		



***** histograms in percent
foreach x in Furnace AirSeal Attic WallIns Window {
		count if Real_tot_act`x'!=. & conthomes==1
		sca totsamp = r(N)
		gen pct_`x'_conthomes = 100*freq_bin`x'_conthomes/totsamp if conthomes==1
}


	count if conthomes==1
	sca totsamp = r(N)
	foreach x in WallIns {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & conthomes==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_conthomes newbin`x' if newbin`x'!=.  & conthomes==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(150) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)80, format(%8.0fc)) yscale(range(0 80)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///								
								xscale(range(-300 2100)) xlabel(-300(300)2100, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_conthomes.png", replace width(2500)
	}

	
	
	foreach x in Furnace {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & conthomes==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_conthomes newbin`x' if newbin`x'!=.  & conthomes==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(150) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-300 3000)) xlabel(-300(300)3000, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_conthomes.png", replace width(2500)
	}

	foreach x in AirSeal {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & conthomes==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_conthomes newbin`x' if newbin`x'!=.  & conthomes==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(50) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-100 900)) xlabel(-100(100)900, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_conthomes.png", replace width(2500)
	}

	foreach x in Attic {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & conthomes==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_conthomes newbin`x' if newbin`x'!=.  & conthomes==1,  ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(100) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-300 2400)) xlabel(-300(300)2400, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_conthomes.png", replace width(2500)
	}

	foreach x in Window {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & conthomes==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_conthomes newbin`x' if newbin`x'!=.  & conthomes==1,  ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(100) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab's") ///
								xscale(range(-200 2000)) xlabel(-200(200)2000, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_conthomes.png", replace width(2500)
	}



********************************************************************************
**** histograms of spending for best versus worst contractors

forval i = 1/5 {
	gen samp_contquint`i' = 1 if contqual_quintile==`i' & conthomes==1
	bysort samp_contquint`i' Household : gen uniq_contquint`i' = 1 if _n==1
	replace uniq_contquint`i' = . if samp_contquint`i'==.
}


forval i = 1/5 {
	foreach x of varlist binFurnace binAirSeal binAttic binWallIns binWindow {
	bysort new`x' uniq_contquint`i' : gen count`x' = _n if uniq_contquint`i'==1
	bysort new`x' : egen freq_`x'_contquint`i' = max(count`x') if uniq_contquint`i'==1
	drop count`x'
	}		
}

***** histograms in percent
forval i = 1/5 {
	foreach x in Furnace AirSeal Attic WallIns Window {
		count if Real_tot_act`x'!=. & uniq_contquint`i'==1
		sca totsamp = r(N)
		gen pct_`x'_contquint`i' = 100*freq_bin`x'_contquint`i'/totsamp if uniq_contquint`i'==1
	}
}


forval i = 1/5 {
	count if uniq_contquint`i'==1
	sca totsamp = r(N)
	foreach x in WallIns {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & uniq_contquint`i'==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_contquint`i' newbin`x' if newbin`x'!=.  & uniq_contquint`i'==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(150) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)80, format(%8.0fc)) yscale(range(0 80)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-300 2100)) xlabel(-300(300)2100, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_contquint`i'.png", replace width(2500)
	}

	foreach x in Furnace {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & uniq_contquint`i'==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_contquint`i' newbin`x' if newbin`x'!=.  & uniq_contquint`i'==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(150) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-300 3000)) xlabel(-300(300)3000, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_contquint`i'.png", replace width(2500)
	}

	foreach x in AirSeal {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & uniq_contquint`i'==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_contquint`i' newbin`x' if newbin`x'!=.  & uniq_contquint`i'==1, ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(50) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-100 900)) xlabel(-100(100)900, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_contquint`i'.png", replace width(2500)
	}

	foreach x in Attic {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & uniq_contquint`i'==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_contquint`i' newbin`x' if newbin`x'!=.  & uniq_contquint`i'==1,  ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(100) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab'") ///
								xscale(range(-300 2400)) xlabel(-300(300)2400, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_contquint`i'.png", replace width(2500)
	}

	foreach x in Window {
	sort newbin`x'
	local varlab : variable label Real_tot_act`x'
	count if Real_tot_act`x'!=. & uniq_contquint`i'==1
	sca `x'_samp = r(N)
	sca `x'_percsamp = round((`x'_samp/totsamp)*100, 0.1)
	local samppercent: display %2.1f `x'_percsamp
	twoway bar pct_`x'_contquint`i' newbin`x' if newbin`x'!=.  & uniq_contquint`i'==1,  ///
								fcolor(blue%30) fintensity(inten30) lwidth(none) barwidth(100) ///
								graphregion(color(white)) bgcolor(white) ///
								bcolor(blue%30) ///
								ylabel(0(10)40, format(%8.0fc)) yscale(range(0 40)) ///
								ytitle("Percent of Homes") xtitle("Amount Spent (US$) on `varlab's") ///
								xscale(range(-200 2000)) xlabel(-200(200)2000, valuelabels labsize(small) angle(45) format(%8.0fc))
	graph export "$maindir/Results/Figures/savpcthist_`x'_contquint`i'.png", replace width(2500)
	}
}





********************************************************************************	
********************************************************************************
***** simulations - how the contractor quality differences affect the wedge

**** quality difference between contractors
sum contqual_pred_new if contsamp==1, detail
sca best5 = r(p5)
sca best10 = r(p10)
sca best25 = r(p25)
sca best50 = r(p50)
sca worst75 = r(p75)
sca worst90 = r(p90)
sca worst95 = r(p95)

foreach x in 5 10 25 50 {
	gen all_tobest`x' = best`x' - contqual_pred_new if contqual_pred_new!=.
}

***** simulations - how the contractor quality differences affect the wedge

foreach x in 5 10 25 50 {
	gen gap_pred_all`x' = savings_gap_pct
	replace gap_pred_all`x' = savings_gap_pct + all_tobest`x' if all_tobest`x'!=.
}


*** finally comparing savings predictions
sum savings_gap_pct if contsamp==1
sca meansave = r(mean)
sca observations = r(N)
matrix contpred_mat = (meansave, 0)

foreach x in 50 25 10 5 {
	sum gap_pred_all`x' if contsamp==1
	sca mean_all`x' = `r(mean)'
	sca pct_all`x' = 100*(mean_all`x'-meansave)/meansave

	matrix contpred_mat`x' = (mean_all`x' , pct_all`x') 
}

drop all_tobest5-gap_pred_all50

***** bootstrapping simulation for effect of contractor quality

matrix bootscontpred_mat = (.,.)
matrix bootscontpred_mat50 = (.,.)
matrix bootscontpred_mat25 = (.,.)
matrix bootscontpred_mat10 = (.,.)
matrix bootscontpred_mat5 = (.,.)

forval i = 1/200 {

	di "starting iteration `i'"

	**** quality difference between contractors
	qui: sum contqualpred_boots`i' [fweight = weight_boots`i'], detail
	sca best5 = r(p5)
	sca best10 = r(p10)
	sca best25 = r(p25)
	sca best50 = r(p50)

	foreach x in 5 10 25 50 {
		gen all_tobest`x' = best`x' - contqualpred_boots`i' if contqualpred_boots`i'!=.
	}

	***** simulations - how the contractor quality differences affect the wedge

	foreach x in 5 10 25 50 {
		gen gap_pred_all`x' = savings_gap_pct`i'
		replace gap_pred_all`x' = savings_gap_pct`i' + all_tobest`x' if all_tobest`x'!=.
	}


	*** finally comparing savings predictions
	qui: sum savings_gap_pct`i' [fweight = weight_boots`i'] if contqualpred_boots`i'!=.
	sca meansave = r(mean)
	matrix tempmat = (meansave, 0)
	matrix bootscontpred_mat = (bootscontpred_mat \ tempmat)

	foreach x in 50 25 10 5 {
		qui: sum gap_pred_all`x' [fweight = weight_boots`i'] if contqualpred_boots`i'!=.
		sca mean_all`x' = `r(mean)'
		sca pct_all`x' = 100*(mean_all`x'-meansave)/meansave
		
		matrix tempmat1 = (mean_all`x' , pct_all`x')
		matrix bootscontpred_mat`x' = (bootscontpred_mat`x' \ tempmat1)
	}
	
	drop all_tobest5-gap_pred_all50

}

matrix contpred_mat_var = bootscontpred_mat[2..., 1...]
mata : st_matrix("contpred_mat_var", mm_colvar(st_matrix("contpred_mat_var")))
matrix contpred_mat_var = diag(contpred_mat_var)
matrix conpredcoefs = contpred_mat
matrix colnames conpredcoefs = "pp" "percent"
matrix colnames contpred_mat_var = "pp" "percent"
matrix rownames contpred_mat_var = "pp" "percent"

bootstrapres conpredcoefs contpred_mat_var observations
est sto baselinecont

foreach x in 50 25 10 5 {
	matrix contpred_mat_var = bootscontpred_mat`x'[2..., 1...]
	mata : st_matrix("contpred_mat_var", mm_colvar(st_matrix("contpred_mat_var")))
	matrix contpred_mat_var = diag(contpred_mat_var)
	matrix conpredcoefs = contpred_mat`x'
	matrix colnames conpredcoefs = "pp" "percent"
	matrix colnames contpred_mat_var = "pp" "percent"
	matrix rownames contpred_mat_var = "pp" "percent"
	bootstrapres conpredcoefs contpred_mat_var observations
	est sto contres`x'
}


esttab baselinecont contres50 contres25 contres10 contres5 ///
				using "$maindir/Results/Tables/simulations_contractors.tex", ///
				replace b(3) se(3) nostar nonotes nonumbers ///
				mlabels("Baseline" "50%" "25%" "10%" "5%" ) ///
				stats(N, fmt(%9.0gc) labels("Observations"))

				
				
				
**********************************************************				
******* model considering interactions between contractor quality and spending in measures

* omitted category will be 20% - 80%
gen oldcontqual_quintile = contqual_quintile
replace contqual_quintile = 3 if contqual_quintile==2 | contqual_quintile==4

reg savings_gap_pct ///
		ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID ib3000.BlowerPreBins ib1.builddate ///
		i.nstories ib1500.binSqft ///
		ib3.contqual_quintile##ib1.binWallIns ///
		ib3.contqual_quintile##ib1.binAirSeal ///
		ib3.contqual_quintile##ib1.binFurnace ///
		ib3.contqual_quintile##ib1.binAttic ///
		ib3.contqual_quintile##ib1.binWindow ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & contsamp==1
		
		
est sto decompose_gap_contquality
gen decompsamp = e(sample)
matrix decomp_coefs = e(b)
sca nobs = e(N)

drop contqual_quintile
rename oldcontqual_quintile contqual_quintile


*** bootstrapping regression to decompose gap interacted with quality
matrix bootsmat = decomp_coefs
forval j = 1/200 {
	matrix tempmat1 = decomp_coefs
	matrix tempmat2 = J(1, colsof(decomp_coefs), .)
	matrix tempmat = (tempmat1 \ tempmat2)
	
	di "starting iteration `j'"
	
	qui : bysort arch_contID ProgramYear : egen boots_contqual = mean(contqualpred_boots`j')

	** different ranking by PY
	forval i = 2010/2015 {
		xtile bootsquint_py`i' = boots_contqual [fweight = weight_boots`j'] if decompsamp==1 & ProgramYear==`i', n(5)
	}

	qui : gen boostrapquint = bootsquint_py2010 if ProgramYear==2010
	drop bootsquint_py2010
	forval i = 2011/2015 {
		qui : replace boostrapquint = bootsquint_py`i' if ProgramYear==`i'
		drop bootsquint_py`i'
	}
	qui : replace boostrapquint = 3 if boostrapquint==2 | boostrapquint==4
	rename contqual_quintile oldcontqual_quintile
	rename boostrapquint contqual_quintile
	
	qui : reg savings_gap_pct`j' ///
		ib1.binFurnace ib1.binAttic ib1.binWallIns ib1.binAirSeal ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib70.binheatkbtu i.heatworks ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		i.white i.black i.haselderly i.hasminor ///
		ib1.roundAtticR  i.ProgramYear ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		i.regionID ib3000.BlowerPreBins ib1.builddate ///
		i.nstories ib1500.binSqft ///
		ib3.contqual_quintile##ib1.binWallIns ///
		ib3.contqual_quintile##ib1.binAirSeal ///
		ib3.contqual_quintile##ib1.binFurnace ///
		ib3.contqual_quintile##ib1.binAttic ///
		ib3.contqual_quintile##ib1.binWindow ///
		[fweight = weight_boots`j'] if decompsamp==1
	
	matrix coefmat = e(b)
	
	local colnames : colnames bootsmat
	foreach x of local colnames {
		local wherecol2 = colnumb(coefmat, "`x'")
		local wherecol1 = colnumb(tempmat, "`x'")
		matrix tempmat[2, `wherecol1'] =  coefmat[1, `wherecol2']
	}	
	

	matrix tempmat = tempmat[2..., 1...]
	matrix bootsmat = (bootsmat \ tempmat)	
		
	drop boots_contqual
	drop contqual_quintile
	rename oldcontqual_quintile contqual_quintile
}
matrix bootsmat = bootsmat[2..., 1...]


mata : st_matrix("gapinter_cont_var", mm_colvar(st_matrix("bootsmat")))
matrix gapinter_cont_var = diag(gapinter_cont_var)

matrix gapinter_cont_coefs = decomp_coefs
local colnames : colnames gapinter_cont_coefs
matrix colnames gapinter_cont_var = `colnames'
matrix rownames gapinter_cont_var = `colnames'

bootstrapres gapinter_cont_coefs gapinter_cont_var nobs
est sto decompose_gap_contquality
parmest, saving("$maindir/Results/Model_Outputs/decompose_contquality_inter.dta", replace)
				

***********************************************************				
** simulations - how the contractor quality differences affect the wedge, eliminating the interactions effects for Wall Insulation

**** quality difference between contractors
qui: sum contqual_pred_new if decompsamp==1 , detail
sca best5 = r(p5)
sca best10 = r(p10)
sca best25 = r(p25)
sca best50 = r(p50)


foreach x in 5 10 25 50 {
	qui : gen all_tobest`x' = best`x' - contqual_pred_new if contqual_pred_new!=.
}

***** simulations - how the contractor quality differences affect the wedge
foreach x in 5 10 25 50 {
	qui : gen gap_pred_all`x' = savings_gap_pct
	qui : replace gap_pred_all`x' = savings_gap_pct + all_tobest`x' if all_tobest`x'!=.
}
	
est restore decompose_gap_contquality	
matrix coefs = e(b)
forval i = 300(300)2400 {
	local wherecol = colnumb(coefs,"5.contqual_quintile#`i'.binWallIns")
	sca inter`i' = coefs[1,`wherecol']
}

foreach x in 5 10 25 50 {
	forval i = 300(300)2400 {
		qui : replace gap_pred_all`x' = savings_gap_pct + all_tobest`x' - inter`i' if all_tobest`x'!=. & contqual_quintile==5 & binWallIns==`i'
	}
}


*** finally comparing savings predictions
qui: sum savings_gap_pct if decompsamp==1
sca meansave = r(mean)
matrix contpred_mat = (meansave, 0)

foreach x in 50 25 10 5 {
	sum gap_pred_all`x' if decompsamp==1
	sca mean_all`x' = `r(mean)'
	sca pct_all`x' = 100*(mean_all`x'-meansave)/meansave
	
	matrix contpred_mat`x' = (mean_all`x' , pct_all`x') 
}
drop all_tobest5-gap_pred_all50	
	

	
************* bootstrapping 
matrix bootscontpred_mat = (.,.)
matrix bootscontpred_mat50 = (.,.)
matrix bootscontpred_mat25 = (.,.)
matrix bootscontpred_mat10 = (.,.)
matrix bootscontpred_mat5 = (.,.)

forval j = 1/200 {

	di "starting iteration `j'"
	
	qui : bysort arch_contID ProgramYear : egen boots_contqual = mean(contqualpred_boots`j')

	** different ranking by PY
	forval i = 2010/2015 {
		xtile bootsquint_py`i' = boots_contqual [fweight = weight_boots`j'] if decompsamp==1 & ProgramYear==`i', n(5)
	}

	qui : gen boostrapquint = bootsquint_py2010 if ProgramYear==2010
	drop bootsquint_py2010
	forval i = 2011/2015 {
		qui : replace boostrapquint = bootsquint_py`i' if ProgramYear==`i'
		drop bootsquint_py`i'
	}
	qui : replace boostrapquint = 3 if boostrapquint==2 | boostrapquint==4
	rename contqual_quintile oldcontqual_quintile
	rename boostrapquint contqual_quintile
	
	**** quality difference between contractors
	qui: sum contqualpred_boots`j' [fweight = weight_boots`j'], detail
	sca best5 = r(p5)
	sca best10 = r(p10)
	sca best25 = r(p25)
	sca best50 = r(p50)

	foreach x in 5 10 25 50 {
		qui : gen all_tobest`x' = best`x' - contqualpred_boots`j' if contqualpred_boots`j'!=.
	}

	***** simulations - how the contractor quality differences affect the wedge
	
	foreach x in 5 10 25 50 {
		qui : gen gap_pred_all`x' = savings_gap_pct`j'
		qui : replace gap_pred_all`x' = savings_gap_pct`j' + all_tobest`x' if all_tobest`x'!=.
	}
	
	
	matrix coefs = bootsmat[`j', 1...]
	forval i = 300(300)2400 {
		local wherecol = colnumb(coefs,"5.contqual_quintile#`i'.binWallIns")
		sca inter`i' = coefs[1,`wherecol']
	}

	foreach x in 5 10 25 50 {
		forval i = 300(300)2400 {
			qui : replace gap_pred_all`x' = savings_gap_pct`j' + all_tobest`x' - inter`i' if all_tobest`x'!=. & contqual_quintile==5 & binWallIns==`i'
		}
	}


	*** finally comparing savings predictions
	qui: sum savings_gap_pct`j' [fweight = weight_boots`j'] if contqualpred_boots`j'!=.
	sca meansave = r(mean)
	matrix tempmat = (meansave, 0)
	matrix bootscontpred_mat = (bootscontpred_mat \ tempmat)

	foreach x in 50 25 10 5 {
	qui:	sum gap_pred_all`x' [fweight = weight_boots`j'] if contqualpred_boots`j'!=.
		sca mean_all`x' = `r(mean)'
		sca pct_all`x' = 100*(mean_all`x'-meansave)/meansave

		matrix tempmat1 = (mean_all`x' , pct_all`x')
		matrix bootscontpred_mat`x' = (bootscontpred_mat`x' \ tempmat1)
	}
	
	drop all_tobest5-gap_pred_all50
	drop boots_contqual
	drop contqual_quintile
	rename oldcontqual_quintile contqual_quintile
}

matrix contpred_mat_var = bootscontpred_mat[2..., 1...]
mata : st_matrix("contpred_mat_var", mm_colvar(st_matrix("contpred_mat_var")))
matrix contpred_mat_var = diag(contpred_mat_var)
matrix conpredcoefs = contpred_mat
matrix colnames conpredcoefs = "pp" "percent"
matrix colnames contpred_mat_var = "pp" "percent"
matrix rownames contpred_mat_var = "pp" "percent"
bootstrapres conpredcoefs contpred_mat_var observations
est sto baselinecont

foreach x in 50 25 10 5 {
	matrix contpred_mat_var = bootscontpred_mat`x'[2..., 1...]
	mata : st_matrix("contpred_mat_var", mm_colvar(st_matrix("contpred_mat_var")))
	matrix contpred_mat_var = diag(contpred_mat_var)
	matrix conpredcoefs = contpred_mat`x'
	matrix colnames conpredcoefs = "pp" "percent"
	matrix colnames contpred_mat_var = "pp" "percent"
	matrix rownames contpred_mat_var = "pp" "percent"
	bootstrapres conpredcoefs contpred_mat_var observations
	est sto contres`x'
}


esttab baselinecont contres50 contres25 contres10 contres5 ///
				using "$maindir/Results/Tables/simulations_contractors_dropinter.tex", ///
				replace b(3) se(3) nostar nonotes nonumbers ///
				mlabels("Baseline" "50%" "25%" "10%" "5%" ) ///
				stats(N, fmt(%9.0gc) labels("Observations"))
				






********************************************************************************
* graph correlating contractor quality and the wedge by measures

*** for wall insulation
clear all
use "$maindir/Results/Model_Outputs/decompose_contquality_inter.dta"


keep if regexm(parm,"WallIns")==1
drop if regexm(parm, "coarsebinWallIns")==1


gen contqualbin = ""
replace contqualbin = "1" if regexm(parm, "1o.cont")==1 | regexm(parm, "1b.cont")==1 | regexm(parm, "1.cont")==1 
replace contqualbin = "2" if regexm(parm, "2o.cont")==1 | regexm(parm, "2b.cont")==1 | regexm(parm, "2.cont")==1 
replace contqualbin = "3" if regexm(parm, "3o.cont")==1 | regexm(parm, "3b.cont")==1 | regexm(parm, "3.cont")==1 
replace contqualbin = "4" if regexm(parm, "4o.cont")==1 | regexm(parm, "4b.cont")==1 | regexm(parm, "4.cont")==1 
replace contqualbin = "5" if regexm(parm, "5o.cont")==1 | regexm(parm, "5b.cont")==1 | regexm(parm, "5.cont")==1 
destring contqualbin, replace


gen measurebin = parm
replace measurebin = subinstr(measurebin, "1.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "1o.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "2.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "2o.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "3.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "3b.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "4.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "4o.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "5.contqual_quintile#", "", .)
replace measurebin = subinstr(measurebin, "5o.contqual_quintile#", "", .)

replace measurebin = subinstr(measurebin, ".binWallIns", "", .)
replace measurebin = subinstr(measurebin, "b", "", .)
replace measurebin = subinstr(measurebin, "o", "", .)

destring measurebin, replace

drop if regexm(parm, "3b.contqual")==1
replace contqualbin = 3 if contqualbin==.


sort contqualbin measurebin
by contqualbin : gen id = _n



** Holm correction for multiple hypothesis testing
sort p
multproc, puncor(0.05) pvalue(p) method(holm) critical(newp_holm) reject(reject_holm)
multproc, puncor(0.05) pvalue(p) method(bonferroni) critical(newp_bonfer) reject(reject_bonfer)
multproc, puncor(0.05) pvalue(p) method(hochberg) critical(newp_hoch) reject(reject_hoch)
multproc, puncor(0.05) pvalue(p) method(krieger) critical(newp_krieg) reject(reject_krieg)
multproc, puncor(0.05) pvalue(p) method(liu2) critical(newp_liu) reject(reject_liu)
multproc, puncor(0.05) pvalue(p) method(simes) critical(newp_simes) reject(reject_simes)

replace reject_simes = 0 if p==.


forval i = 1(2)5{
eclplot estimate min95 max95 id if contqualbin==`i', ///
    supby(reject_simes) ///
	estopts1(mcolor(blue) msymbol(Th))  ///
	estopts2(mcolor(red) msymbol(S)) ///
	ciopts1(lcolor(blue)) ///
	ciopts2(lcolor(red)) yline(0) ///
	legend(off) ///
	xtitle("Amount Spent (US$) on Wall Insulation") ///
	graphregion(color(white)) bgcolor(white) ///
	ytitle("Performance Wedge (%)" "(compared to zero spending)") yline(0) ///
	xlabel(1 "0" 2 "1-300" 3 "301-600" 4 "601-900" 5 "901-1200" 6 "1201-1500" ///
	7 "1501-1800" 8 "1801-2100" 9 ">2100" , ///
	labsize(small) angle(45))

graph export "$maindir/Results/Figures/wedge_wallins_contqual`i'.png", replace width(2500)
}
